A continuous-time solver for quantum impurity models 
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We present a new continuous time solver for quantum impurity models such as those relevant 
to dynamical mean field theory. It is based on a stochastic sampling of a perturbation expansion 
in the impurity-bath hybridization parameter. Comparisons to quantum Monte Carlo and exact 
diagonalization calculations confirm the accuracy of the new approach, which allows very efficient 
simulations even at low temperatures and for strong interactions. As examples of the power of the 
method we present results for the temperature dependence of the kinetic energy and the free energy, 
enabling an accurate location of the temperature-driven metal-insulator transition. 

PACS numbers: 71.10.-w, 71.10.Fd, 71.28.+d, 71.30.+h 
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Numerical computation of dynamical properties of 
strongly correlated fermion systems is one of the fun- 
damental challenges in condensed matter physics. The 
exponential growth of the Hilbert space renders exact di- 
agonalizations impossible except for very small systems 
When applied directly to lattice models, auxiliary 
field Monte Carlo methods encounter severe difficulties 
with the fermionic sign problem 0. The density matrix 
reormalization group 3] is powerful for extracting ground 
state properties of one dimensional systems, but exten- 
sions to higher dimensions and dynamics have proven 
difficult. Over the last decade, a promising new "dynam- 
ical mean field" (DMFT) approach has been developed 
As shown by Georges, Kotliar and other workers 
if the momentum dependence of the self-energy 
is neglected, E(p, uS) — > £(oj), then the solution of the 
lattice model may be obtained from the solution of a 
quantum impurity model plus a self-consistency condi- 
tion. 

A quantum impurity model represents an atom or 
molecule embedded in a host medium; in addition to 
their relevance to DMFT calculations these models are of 
intrinsic interest and important to nano-science as rep- 
resentations of single-molecule conductors. A quantum 
impurity model consists of a set of levels a (with cre- 
ation operators "01), populated by electrons which inter- 
act via general four-fermion interactions JJ abcd (for exam- 
ple density-density or spin exchange terms). The levels 
a are hybridized to electronic continua ( "bath orbitals" ) 
representing the degrees of freedom of the host material. 
These bath degrees of freedom may be integrated out, 
and the impurity model specified by a partition function 
Z = TrT T e with effective action S = Sf + S\ oc , where 



dTdT^ a {T)F a {T-T')^ a {T') 



(1) 



S'loc — 



dr(e ab ^ b - U abcd ^l^ d ). (2) 



The function F is determined by the hybridization and 



bath density of states. It describes transitions from the 
impurity into the bath and back and is related to the 
mean field function Gq 1 Q by G^ 1 ^) = iu> + /i — 
F{—iuj). A solution of the quantum impurity model 
amounts to evaluating the Green functions G(r) = 
-{T t i/}(t)^(0)) for given F, e and U . The dynamical 
mean field procedure involves a self-consistent determi- 
nation of the hybridization function F, in which the com- 
putation of G is time critical 7] . 

Perhaps the most commonly used technique is the 
Hirsch-Fye method |9j, in which a (discrete) Hubbard- 
Stratonovich transformation is used to decouple the in- 
teraction part, leading to determinants which give the 
weights associated with the configurations of the auxil- 
iary fields, which are then sampled by a Monte Carlo 
procedure. The method requires a discretization of the 
imaginary time interval [0, 0] into N equal slices and the 
evaluation of determinants of N x N matrices. At low 
temperatures and strong correlations, the Green func- 
tions have a highly non-uniform time dependence, drop- 
ping very steeply near the edge points r = 0, 0, and 
varying more slowly in between. The rapid drop means 
that a very large N (~ 5U/3) is required, limiting the 
utility of the algorithm for physically relevant temper- 
atures. Furthermore, useful decouplings of the general 
interactions of multiorbital models are lacking. Another 
approach is an exact diagonalization method [ToL ITlj in 
which the continuum of energies in the impurity model is 
replaced by a discrete set of levels. This method suffers 
from systematic errors associated with the discrete level 
spacing, and becomes impractical in the multiorbital case 
because of the rapid growth of the Hilbert space. 

In this paper we present a new, continuous time ap- 
proach which resolves many of these issues. The idea, 
introduced in the context of bosonic field theories by 
Prokof'ev et al. ' s t° define a formal perturbation 
expansion for the partition function Z . Certain collec- 
tions of terms in this series are then sampled by a Monte 
Carlo procedure. In important prior work, Rubtsov and 
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co-workers ^jj performed such a stochastic evaluation by 
expanding in the interaction term. In this paper we show 
how to formulate a determinental Monte Carlo algorithm 
by expanding in the hybridization term Sp while treat- 
ing the interactions exactly. This strong-coupling ap- 
proach becomes particularly powerful at the strong inter- 
actions characteristic of systems of present day interest 
(such as high-T c cuprates), because the perturbation or- 
der actually decreases with increasing U. Our method al- 
lows access to low enough temperatures that both ground 
state properties and the leading temperature dependen- 
cies may be determined, even at very strong couplings, 
providing new information unavailable by other methods. 

We illustrate the procedure first for noninteracting 
spinlcss fermions, and then proceed to the one-orbital 
Hubbard model, reserving a general discussion for a fu- 
ture publication Q . Because the spinless fermion model 
is diagonal in the occupation number basis, we have for 
Z F = Z(ji = 0, U = 0, cr) 

ZF=2 + J2 dT l dT l dT 2 dT 2--- 

k—1 i i T 2 



drt 



dT e k Z^(rl,Tl;Ti,Ti;...;r^), (3) 



with -^(t^, rf ; t| , r|; . . . ; t|, r k ) denoting the collection 
of k\ diagrams containing k segments of occupied fermion 
states, with starting point r? and end point rf , i — 
1, . . . , k. We view the interval [0, f3] as a circle and take 
care of the trace over occupation number states by includ- 
ing diagrams in which the last fermion line winds around 
the circle. This is denoted by an upper integral bound 
or. An illustration of these diagrams for k = 3 is shown 
in Fig. ^ Each end point is connected to a starting 
point by a dashed line representing F(t^ — r*). The 
collection of diagrams evaluates to 
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.;r fe s ) r fc e )=detM- 1 ^ 1 (4) 
- <), (5) 



with = — 1 if r k < t( (the last segment winds around 
the circle ) and +1 otherwise. One also has to sample 
the "full line" and the "empty line" states, correspond- 
ing to 1 or particle in the whole interval < t < (3. 
We use three types of Monte Carlo updates: (i) insertion 
and removal of a segment, (ii) insertion and removal of 
an anti-segment, and (iii) shift of the segment end point. 
The first two are required for ergodicity and the third 
enhances the sampling efficiency. Combining the fc! di- 
agrams into a determinant is essential. Some diagrams 
have negative weight and a "worm" -type algorithm which 
sampled individual diagrams would run into a bad sign 
problem. Yet we found in all of our computations that 
the determinant remains positive. 
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FIG. 1: Illustration of Z$ ' (Ti,Ti;Ti,T2;T3,T§) in terms of 
diagrams. A full line corresponds to an occupied fermion 
state. Each end point of a segment is connected to a 
starting point by a dashed line representing F(r„ — t%). 
The combined weight of the diagrams can be expressed as a 
determinant. 



Suppose we start with a collection of segments = 
r fc' r fe) an d want to insert a new seg- 
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ment s at a randomly chosen f s . If f s happens to lie 
on a segment of s k , the move to the new configuration 
Sk+i = s k + s is rejected, otherwise it has to satisfy a 
detailed balance condition. If the length I of the new 
segment is chosen randomly in the interval [0, Z max ]j with 
Zmax determined by f s and the collection s k of existing 
segments, and we propose to remove this new segment 
in the new configuration Sk+i with probability l/(fe + l), 
the condition is 



Pins (^r 



Prem (^d 



Z[ +1 (s k+1 ) [31 



max ^l/j, 



1 



(6) 



The result for anti-segments is similar and for shift up- 
dates from a configuration s to a configuration s (chang- 
ing the length of some segment from I to a randomly 
chosen I G [0, l maK ]) one obtains 



p(s -> S) 



zE{s) 



(7) 



A procedure analogous to that in Ref. 13] is used to 
calculate the determinant ratios and the new enlarged 
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FIG. 2: Green functions for n = 1, U/t = 3.5^2, fit = 
400, 200, 31.4 and 20. Lines without symbols (upper and right 
axes) show G(r) on a semi-log scale over the wide time in- 
terval [(3/2, (3] revealing marked differences between metallic 
(pt = 200, 400) and insulating (fit = 20, 31.4) solutions. Lines 
with symbols (lower and left axes) show the same data on a 
linear scale in the very narrow r range [0,/3/2000], revealing 
the accurate representation of the rapid drop of G(r). 



(reduced) matrices in a time 0(k 2 ). We store and ma- 
nipulate M, the inverse of Eq. (JSJ, because M allows easy 
access to the determinant ratios in Eqs. JHJ) and and 
is required for measuring the Green function, since 



k k 

P 8=1 j = l 

t' > 
-5(t -t'-(3) t' < ' 



(9) 



The end points G(0) and G(/3) can be measured accu- 
rately from the average total length of the segments. 

In the form given here, the algorithm generalizes 
straightforwardly to any model with interaction terms 
which are diagonal in an occupation number basis (for 
models with exchange, see Ref. 01)- O ne simply in- 
troduces one collection of segments for each spin/orbital 
state, and the weight of a configuration now also de- 
pends on the segment overlap. For example, in the one- 
orbital Hubbard model with on-site interaction U, there 
is one collection of segments for spin up and one for spin 
down, while in Eqs. © and J7J one has to add a factor 
exp(— S ov U) on the right hand side, where 8 v denotes 
the change in overlap between up and down segments. 

We have used the new method to study the param- 
agnetic phase of the Hubbard model with semicircular 
density of states of bandwidth 4t, for interactions of the 
order of the Mott critical value U C 2 and temperatures 
as low as fit = 400. For this model the self-consistency 
condition reduces to F(t) = t 2 G(— t). Simulations for 
temperatures down to (3t ~ 50 can be run on a laptop. 
For calculations at fit = 400, we typically used 10 CPU 
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FIG. 3: Probability p(k) for a configuration with k segments 
plotted for different interaction strengths U/t for fit = 100 
and half-filling. The peak position shifts to lower values of 
k as U/t is increased. The inset compares the scaling of the 
matrix size with U to Hirsch-Fye (~ 5/3(7) and the method of 
Ref. 0] (« 0.5/3(7). 



hours for each iteration in order to accurately resolve the 
short- and long-time behavior. 

Figure [3 shows the impurity model Green function for 
U/t = 3.5^, (3t = 20, 31.4, 200 and 400 and n = 1 (half 
filling). The lower two temperatures are out of reach of 
the Hirsch-Fye algorithm. We collected the data on a grid 
of 10 4 points for 0t = 200,400 and 10 3 points for f3t = 
20,31.4. The lines with symbols show that the method 
accurately captures the steep short-time drop of G; the 
lines without symbols demonstrate clearly the difference 
in long-time behavior between the insulating (high-T) 
and metallic (low-T) solutions. 

Despite the almost perfect resolution, the typical size, 
k, of the matices, M, which are generated during the 
simulation remains reasonable even at low temperatures. 
This property explains the superior performance of the 
strong-coupling expansion method. Figure 01 shows the 
probability distribution p{k) for fit — 100 and different 
values of the interaction strength. While the peak value 
of the distribution is proportional to (3, it shifts to lower 
order as the interaction strength is increased, in contrast 
to Hirsch-Fye or the method of Ref. [l3j , where the ma- 
trix size scales approximately as 5/3(7 and 0.5/3(7, respec- 
tively. The inset of Fig.|3shows that the linear size of the 
matrix in our method can easily be a factor 100 smaller 
than in a Hirsch-Fye calculation or a factor 10 smaller 
than in the weak-coupling approach of Ref. |l3| . The cu- 
bic scaling of the computational effort with matrix size 
implies a dramatically improved efficiency at couplings 
of the order of the Mott critical value, making low T 
behavior accessible. 

To verify the accuracy of the method we show in Fig. 01 
the kinetic energy K — 2t 2 J (ItG{t)G{—t) obtained 
via the new approach, the exact diagonalization method 
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FIG. 4: Kinetic energy obtained using the indicated impurity 
solvers plotted as a function of temperature for U/t = 4 and 
U/t — 3.5\/2 ~ 4.95. The former value corresponds to a 
metallic phase, while for U/t — 3.5V2, the system undergoes 
a metal-insulator transition. 
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FIG. 5: Points: Kinetic energy as a function of temperature 
for U/t = 5.3. Solid line with error estimates as dotted lines: 
free energy of the metallic state. Dashed line: free energy of 
the insulating state. Inset: quadratic T-dependence of the 
total energy. 



|15| . and the Hirsch-Fye method. The results are plot- 
ted against T 2 to show that the theoretically expected 
Fermi liquid result can be obtained even when the char- 
acteristic energy scales are very low. The new algorithm 
displays the metal-insulator transition clearly, agrees per- 
fectly with exact diagonalization results, and is consistent 
with the Hirsch-Fye data at elevated temperatures. 

Figure 0] shows a discontinuity arising from the in- 
stability of the metallic solution as the temperature 
is increased above a certain value (T/t w 0.026 at 
U/t = 4.95). The instability implies the existence 
of a first order metal-insulator transition, but strictly 
speaking it is a spinodal point. To further demon- 
strate the power of the method we locate the true phase 
boundary by constructing the free energy of the metal- 
lic state. We compute the total internal energy per 
site E(T) — K(T) + U(n^ni), from which we obtain 
the specific heat C(T) = dE(T)/dT and hence the en- 
tropy S(T) = / r dT'(C(T')/T'). We compare this to 
the free energy of the paramagnetic insulating state for 
which E(T) has negligible T-dependence because of the 
gap, but the entropic contribution is —Tin 2. Results 
are illustrated in Fig. for U/t = 5.3, where the in- 
teraction is so large that exact diagonalization data in 
the relevant T regime are not available. The jump in 
the kinetic energy is clearly seen to be a spinodal ef- 
fect, occuring at a higher temperature (T « 0.01 It) 
than the T c — 0.00685(4)i at which the free energies 
cross. The latent heat of this metal-insulator transition 
is T c (S mct (T c ) - S ins ) = 0.00250(5)i. 

In conclusion, we have developed a continuous-time im- 
purity solver, based on a stochastic evaluation of an ex- 
pansion in the hybridization. The new approach, which 
(even away from half-filling) does not appear to suffer 
from a sign problem, is much more efficient at intermedi- 



ate to strong couplings than other available methods. It 
allows to simulate the Hubbard model at previously in- 
accessible temperatures, permitting direct access to the 
asymptotic low energy physics, for example allowing the 
computation of the latent heat across the metal-insulator 
transition. The algorithm opens up for systematic study 
the low temperature properties of strongly interacting 
quantum models relevant to transition metal oxides, ac- 
tinides and other correlated electron materials. 
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